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Abstract 



In this paper we propose a monotone numerical scheme for fully 
nonlinear parabolic PDEs, which includes the quasi-linear PDE asso- 
ciated with a coupled FBSDE as a special case. Our paper is strongly 
. motivated by the remarkable work Fahim, Touzi and Warin [15], and 

stays in the paradigm of monotone schemes initiated by Barles and 
Souganidis [3]. Our scheme weakens a critical constraint imposed by 
[15], especially when the generator of the PDE depends only on the 
CN ! diagonal terms of the hessian matrix. Several numerical examples, up 

to dimension 12, are reported. 
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1 Introduction 



In this paper we are interested in feasible numerical schemes for the following 
fully nonlinear parabolic PDE, especially in high dimensional case: 

- d t u - G(t, x, u, Du, D 2 u) = on [0, T) x R d ; u(T, ■) = g{-) on R d . (1.1) 
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The standard numerical schemes in the PDE literature, e.g. finite difference 
method and finite elements method, work only for low dimensional problems, 
typically d < 3, due to the well known curse of dimensionality. However, in 
many applications, especially in finance, the dimension d can be higher. We 
thus turn to probabilistic approach which is less sensitive to the dimension. 

In the semilinear case, the PDE (1.1) is associated to a Markovian Back- 
ward SDE due to the nonlinear Feynman-Kac formula introduced by Pardoux 
and Peng [23]. Based on the regularity results of BSDEs established by Zhang 
[28], [28] and Bouchard and Touzi [8] proposed the so called Backward Eu- 
ler Scheme for such BSDEs and hence for the associated semilinear PDEs, 
and obtained the rate of convergence. This scheme approximates the BSDE 
by a sequence of conditional expectations, and several efficient numerical 
algorithms have been proposed to compute these conditional expectations, 
notably: Bouchard and Touzi [8], Gobet-Lemor-Waxin [17], Bally-Pages- 
Printems [1]), Bender-Denk [5], Crisan-Manolarakis [12]. There have been 
numerous publications on the subject and the schemes have been extended to 
more general BSDEs, e.g. reflected BSDEs which is appropriate for pricing 
and hedging American options. Typically these algorithms work for 10 or 
higher dimensional problems. 

We intend to numerically solve PDE (1.1) in fully nonlinear case, in par- 
ticular the Hamilton- Jacobi-Bellman equations and the Bellman-Isaacs equa- 
tions which are widely used in stochastic control and in stochastic differential 
games. We remark that this is actually one main motivation of the develop- 
ments of second order BSDEs by Cheridito, Soner, Touzi and Victoir [9], and 
Soner, Touzi and Zhang [25]. Our scheme is strongly inspired by the work 
Fahim, Touzi and Warin [15]. Based on the monotone scheme of Barles and 
Souganidis [3], [15] extended the backward Euler scheme to fully nonlinear 
PDE (1.1). In the case G is convex in (u, Du, D 2 u), they obtained the rate of 
convergence by using the techniques in Krylov [18] and Barles and Jakobson 
[2]. They applied the linear regression method, see e.g. [17], to compute 
the involved conditional expectations, and presented some numerical exam- 
ples up to dimension 5. We remark that the rate of convergence has been 
improved very recently by Tan [26], by using purely probabilistic arguments. 

There is one critical constraint in [15] though. In order to ensure the 
monotonicity of the backward Euler scheme, they assume the lower and up- 
per bounds of G? 7 , the derivative of G with respect to D 2 u, satisfies certain 
constraint. However, when the dimension is high, this constraint implies that 
G 7 is essentially a constant and thus the PDE (1.1) is essentially semilinear, 
see (2.8) for more details. This is of course not desirable in practice. 

The main contribution of this paper is to relax the above constraint. 
In [15] and most papers in the literature of numerical BSDEs, the involved 
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conditional expectations are expressed in terms of Brownian motion. Our 
first simple but important observation is that we may replace the unbounded 
Brownian motion with bounded trinomial tree, which helps to maintain the 
monotonicity of the scheme. We next modify the scheme further, but still in 
the paradigm of monotone scheme, so as to relax the constraint. In the special 
case where G 1 is diagonal, namely G involves D 2 u only through its diagonal 
terms, the above constraint is removed completely. Rate of convergence is 
also obtained. Several numerical examples are presented. In low dimensional 
case, our scheme is comparable to finite difference method and is superior to 
the simulation methods. When G 7 is diagonal, our scheme works well for 12 
dimensional problems. 

We note that PDE (1.1) covers the quasilinear PDEs as a special case, 
which corresponds to a coupled forward backward SDE due to the four step 
scheme of Ma, Protter and Yong [19]. There are only a few papers on nu- 
merical methods for FBSDEs, e.g. Douglas-Ma-Protter [19], Makarov [21], 
Cvitanic and Zhang [11], Delarue and Menozzi [13], Milstein-Tretyakov [22], 
Bender and Zhang [7], and Ma, Shen and Zhao [20]. Most of them deal with 
low dimensional FBSDEs only, except that [7] reported a 10- dimensional nu- 
merical example. However, [7] proved the rate of convergence only for time 
discretization, and the convergence of the linear regression approximation 
is not analyzed theoretically. Our scheme works for FBSDEs as well, espe- 
cially when the diffusion coefficient a is diagonal. A numerical example for 
a 12-dimensional coupled FBSDE is reported. 

We have also presented a few numerical examples which violate our as- 
sumptions and thus the scheme may not be monotone. Numerical results 
show that our scheme still converges. In particular, we note that our current 
theoretical result does not cover the G-expectation, a nonlinear expectation 
introduced by Peng [24]. We nevertheless implement our algorithm to ap- 
proximate a 10-dimensional HJB equation, which includes the G-expectation 
as a special case, and it indeed converges to the true solution. It will be very 
interesting to investigate the convergence of our scheme, or its variations if 
necessary, when the monotonicity condition is violated. We shall leave this 
for our future research. 

The rest of the paper is organized as follows. In section 2 we present some 
preliminaries. In Section 3 we propose our scheme and prove the main con- 
vergence results. Section 4 is devoted to the study of quasilinear PDEs and 
the associated coupled FBSDEs. In Section 5 we discuss how to approximate 
the involved conditional expectations. Finally we present several numerical 
examples in Section 6, up to dimension 12. 
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2 Preliminaries 



Let T > be the terminal time, d > 1 the dimension of the state variable x, 
E> d the set of d x d symmetric matrices. For x,x E~R d and 7, 7 G §> d , denote 

d 

x ■ x := XjXj, |x| := a/x • x, and 7:7:= tr( 77 ), I7I := ^7 : 7. 

i=l 

For x G R d , x T denotes its transpose, and thus xx T G E> d . For 7, 7 G E> d , we 
say 7 < 7 if 7 — 7 is nonnegative definite. For any 7 = [ 7 y] G S d , denote 

-D[ 7 ] := the diagonal matrix whose (i,i)-th component is 7 &. (2.1) 

It is clear that, for any 7, 7 G §> d , 

D[ 7 ]:7 = D[ 7 ]:D[7]=7:D[7]. (2.2) 

Moreover, we use the same notation to denote the zeroes in R d and E> d . 

Our objective is the PDE (1.1), where G : (t,x,y,z,i) G [0,T] x R d x 
R x R d x § d — ^ R and g : x G R d — > R. We shall always assume the following 
standing assumptions: 

Assumption 2.1. (i) G(-, 0,0,0) and g are bounded. 

(ii) G is continuous in t, uniformly Lipschitz continuous in (x, y,z,j), 
and g is uniformly Lipschitz continuous in x 

(Hi) PDE (1.1) is parabolic, that is, G is nondecreasing in 7. 

(iv) The PDE (1.1) has a unique bounded viscosity solution u, and the 
comparison principle for its bounded viscosity solutions holds true. 

For notational simplicity, throughout the paper we assume further that 

G is different iable in (y, z, 7) so that we can use the notations G 7 etc. 

However, we emphasize that all the results in the paper do not rely on this 
additional assumption. For the theory of viscosity solutions, we refer to the 
classical references Crandall, Ishii and Lions [10] and Fleming and Soner [16]. 
Our goal of the paper is to numerically compute the viscosity solution u. In 
their seminal work Barles and Souganidis [3] proposed a monotone scheme 
in an abstract way and proved its convergence by using the viscosity solution 
approach. To be precise, for any t G [0, T) and h > with t + h < T, let 
be an operator on the set of measurable functions : R d ->• R. Now n > 1, 
denote h := — , U := ih, i = 0, 1, • • • , n, and define: 

u h (t n ,x) :=g{x), u h (t, ■) := Tl^tluhiti, •)], te[ti-i,t«), % = n, ■ ■ ■ , 1. (2.3) 

The following convergence result is due to Fahim, Touzi and Warin [15] The- 
orem 3.6, which is based on [3]. 



4 



Theorem 2.2. Let Assumption 2.1 hold. Assume the operator satisfies 
the following four conditions: 

(i) Consistency: for any (t,x) G [0,T)xR d and any ip G C 1,2 ([0, T) xl d ), 

[c + ^(f,x')-T*[[c + ^](f,-)](x') 



lim 

(t',x',/i,c)-»(t,a:,0,0) fl 

= —d t (p(t, x) — G(t, x, if, Dp, D 2 Lp). 

(ii) Monotonicity: T^[</?] < T^[/0] whenever if < if). 

(Hi) Stability: Uh is bounded uniformly in h whenever g is bounded. 

(iv) Boundary condition: for any x G M. d , 

liminf Uh(t',x')= limsup Uh(t',x'). 

(V ,x' ,h)-*(T,x,0) (t' ,x' ,h)->(T,xfl) 

Then Uh converges to u locally uniformly as h — > 0. 

[15] proposed a scheme as follows. Assume there exist a, a G E> d such 
that < |ct 2 < G 7 < \a 2 . Denote 

F(t, x, y, z, 7) := G(t, x, y, z, 7) - ^a 2 : 7, (2.4) 

and define 

where, for a <i-dimensional standard Normal random variable N, 

f>iip(x) := E ip(x + \ZhaN)Ki(N) , i = 0,1,2, (2.6) 

K (iV):=l, K^N) K 2 {N):- 

This scheme satisfies the consistency, and the stability follows from the 
monotonicity. However, to ensure the monotonicity, one needs to assume 
F 1 : g_~ 2 < 1, see the proof of [15] Lemma 3.12. This essentially requires 

[^a 2 - ^q 2 } : a~ 2 < 1, and thus a 2 : a~ 2 < d + 2. (2.7) 

In the case |a 2 = aid, \o 2 = al d , we have 

\<a/a<l + 2- ( 2 - 8 ) 

When d is large, this implies a ~ o. and thus G is essentially semilinear, 
which of course is not desirable in practice. 

Our goal of this paper is to modify the algorithm (2.5)-(2.6) so as to 
relax the above constraint. In particular, in the case that G 7 is diagonal, we 
remove this constraint completely. 
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3 The numerical scheme 

In this section we present our numerical scheme and study its convergence. 
Our scheme involves two parameters < o~$ e § d and < p < 1, which will 
be specified in Section 3.2 below. We shall always denote 

F(t,x,y,z,j) := G(t,x,y,z,-y) - -a^ : 7, G 7 := a 1 G 7 o- _1 . (3.1) 

However, unlike in [15], we emphasize that we do not require \o~q < G 1 . Let 
(Q, J 7 , P) be a probability space, and £ : £7 — > M. d be a random variable such 
that its components £ i( z = 1, • • • , d are independent and have the identical 
distribution: 

p (6 = -^) = |. p te = -^) = !> p(6 = o) = i-p- (3.2) 

This implies that 

E[fc] = E[^ 3 ] = 0, E[tf] = l, E[£ 4 4 ] = i. (3.3) 
We now modify the algorithm (2.5)-(2.6): 

TM(x) := Vfo(x) + hF(t,x,Vfo(x),VMx),Vfo(x)), (3.4) 
where, recalling (2.1), 

Vi V (x) :=E[ V (x + Vha„i)K i (£)\ , i = 0,1,2, (3.5) 
K ® := !,*({) := := '^■-^^^'k 1 . 

One may check straightforwardly that 

^i(£)] = 0, E[K 2 (£)] = 0. (3.6) 
We recall that the approximating solution Uh is defined by (2.3). 

Remark 3.1. If we assume ^<7g < G 7 and set p — | ; i/ien ow scheme is 
obtained by replacing the normal random variable N in (2.6) with trinomial 
random variable £. 27ms m /aci /ias already been mentioned in [15]. 
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3.1 Consistency 

We first justify our scheme by checking its consistency. 

Lemma 3.2. Under Assumption 2.1, T\ satisfies the consistency requirement 
in Theorem 2.2 (i). 

Proof. We first assume (f, x', c) = (t, x, 0) and send h j. 0. Let ip e C 1,2 ([0, T] x 
W d ). Apply the Taylor expansion we have 

ip(t + h, x + Vho~o(,) = (p(t, x) + d t <f(t, x)h + VhDLp(t, x) ■ <7o£ 

+^D 2 ip(t,x) : [a ^[a ^ T + o(h). 

By (3.3) and the independence of one may check straightforwardly that 



V° h cp(t + h,-)(x) = ip(t,x) +d t cp(t,x)h+ \D 2 ^(t,x) : a 2 + o(h); 
V\ip{t + h, -)(x) = Dp(t,x) + o(Vh). 

Moreover, for any i ^ j, 



(3.7) 



E 
E 
E 
E 



0; 



(1 " p)tt T " (1 " 3p)D[tf T ] - 2ph\ = 0; 
6 [(1 - P)Zf - (1 - Z V )DU T \ - 2pl d ] = 
e|[(l ~p)tt T - (1 - 3p)DU T ] - 2pl d ]] = 2(1 -p)S i!i; 
6£i [(1 - P)tt T - (1 - 3p)D[K T ] - 2pl d ]\ = (1 - p)(S id + S_ 



Here Sij is the d x d-matrix whose (i,j)-th component is 1 and all other 
components are 0. Then, denoting A = [a it j] := o- D 2 ip(t, x)a , 



E 



(D 2 y{t,x) : Mhfjao^Kl-p)^ - (1 - 3p)D[^ T ] ~ 2phW 

<j^e[(a : a T ) [(i - p)ee T - (i - zp)Dm T ] - 2 P i d \ 

d 

a, 1 ^ feO [(1 - p)tt T - (1 - 3p)D[K T ] - 2pl d ] 



-1 



*,i=i 

^l-pK-V = 2(1 -p)D 2 <p(t,x), 
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and thus 

V 2 h ip(t + h,-)(x) = D 2 <p(t,x) + o(l). (3.8) 
Combine (3.7) and (3.8), one sees immediately that 

r P(*,3Q-TWV)](s) « ( , v n n2 s 

lim d = —Ot^(t,x) — G(t,x,(f,Dip,Dip). 

h\a h 

The general consistency follows from straightforward extension of the 
above arguments and we omit the details. □ 

3.2 The monotonicity 

To obtain the monotonicity of our scheme, we need to impose the following 
key assumption of the paper. 

Assumption 3.3. There exist < a G S d , scalar constants < a < a and 
< 9 < § such that 

(i) al d < D[GA < al d and D[G 7 ] < (1 + 9)G 7 ; 

(ii) A := a / a < Ag : = sup e i A(p, 60, where 

K 7 ' — 2(1+0) -P- 3 V 7 

A( "- 9): - 1+ 2^(1 + ^-1) - (X9) 

(mj i??/ otherwise reseating o"o (namely multiplying o~q by a scalar con- 
stant) we may set without loss of generality that 

p G \ 2 (i+e) 1 1] suc h ^hat A < X(p, 9) < A 9 ; 

A , _ p(2+30)-9 ^' 10 ' ) 

a ~~ 2 P (A-i)+a p waere a p .— p ( 1+(9 ) ■ 

Remark 3.4. In this remark we provide a few facts concerning our choices 
of parameters which will be used in the proof of next lemma. 

(i) We need p < | so that 1 — 3p, the coefficient of -D[££ T ], is nonnegative. 
For < 9 < |, we have 2 (i+e) < §■ Moreover, for p G [ 2 (i%) > §]> h holds that 
a P > 2p and A(p, 6>) > 1. 

(ii) Note that A is invariant of the scaling of <To and recall that a = a/ A. 
Then our choice of a implies 2a = ^ + (2 — ^f)a. □ 

Lemma 3.5. Let Assumptions 2.1 and 3.3 hold. Then there exists a constant 
h > 0, depending on d, T, 6,p, A, X(p, 9), and the bound and Lipschitz con- 
stants in Assumption 2.1, such that T l h satisfies the monotonicity in Theorem 
2.2 (ii) for all h G (0,h }. 
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Proof. Let <pi < ip 2 be bounded and ip := (p 2 — <fi > 0. Then by (3.4) we 
have, at (t, x), 

T h & 2 ] - T h [(pi] = V° h ^ + h [F y V Q h ^ + F z ■ + F 7 : V\^\ . (3.11) 

Here the terms F y ,F z , F 7 are defined in an obvious way and we emphasize 
that they are deterministic. Plug (3.5) into above equality, then 



= E 
= E 
> E 



i/>{x + VhaoO [I + h[F y + F z ■ K.iO + F, : K 2 (£)] 

I 



ij(x + Vhaot) [l + hF y + VhF z ■ (a ^) 
7p(x + Vho-ot) [1 - CVh + ] , 

thanks to Assumption 2.1, where, 

I := F 1 : (V[(l - p)tt T - (1 - 3p)D[£?] - 2pl d ]a 1 



1 -pL 
(3.12) 



[C 7 - \l d \ : [(1 - p)tt T - (1 - 3p)D[tf T ] - 2pl d ] 
1 . , _ ~ , ~ , 1 



= (1 - P)[G 7 " : (ft ' ) - (1 - 3p)[^[G 7 ] - : (^ ) 
-2ptr(G 7 ) +pd. 

Denote «j := (G 7 )a. Then it follows from Assumption 3.3 (i) that 

i > a-p)[Y^D[G,]-\i4--^e)-(i-mD[G,]-li4--^e 



+ 

-2ptr(G 7 ) +pd 
p(2 + 36)-6 



D[G 7 ];{te)-vW-lMGi)+vd 



pd - p ki + 2 «i - « P «iC 2 



i=l 



Note that a < on < a, a v > 2p by Remark 3.4 (i), and £g takes only values 



and -. Then 

p 



| £v r 1 q/ 

^ + 2a { - a p a£f < (2atA V [- + 2a { ^aA < 2a V - + (2 - -£)a 

p p Vp p 

2A 

= 2a = — — : , 

2p(A- 1) + a p 7 
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thanks to Remark 3.4 (ii) and Assumption 3.3 (iii). Thus 

2A 

I > pd — pd- 



'2p(A- 1) +Op" 

Note that the right side above is decreasing in A. Then, for h small enough, 



i-cVh + - — > i + 



p i — p 



pd — pd- ^ 



2p(X(p,9)-l) + a, 



0. 



where the last equality can be checked straightforwardly by using the defini- 
tion of X(p,8). This, together with (3.12), proves the monotonicity. □ 

Remark 3.6. In this remark we comment on the constraint A < Ag and the 
optimal choice of ao and p. 

(i) We note that A and Ag depend on cr , but not on p. Moreover, they 
are invariant of the scaling of ctq. We shall first choose cq among those with 
unit determinant to maximize and then rescale ao as in Assumption 3.3 
(iii). 

(ii) When 9 = 0, namely G 1 is diagonal, we have Ao = oo. Then we 
remove the constraint (2.8) completely and thus improve the result of [15] 
significantly. In this case we may choose p > small enough so that A < 
\(p, 0) = 1 + ^pzij • We also note that when d = 1 we always have 9 = 0. 

(iii) When (d — 3)9 < 2 and 9^0, the optimal p is 2 -(d-S)o wm ch belongs 

to the interval [ 2 /i + 0\ > §]■ 111 this case Ag = 1 + 8 e(i+e)fi-i) • 

(iv) When (d — 3)9 > 2, the optimal p is |, and in this case Ag = 

j _|_ 3{2-d6) 



2(l+0)(d-l) • 

(v) When 9 > |, Assumption 3.3 is violated. In this case we may always 
set p := | and our algorithm reduces back to [15], by replacing the Brownian 
motion there with trinomial tree. We may easily obtain the bounds (2.7) and 
(2.8) as in [15]. See also Remarks 3.1. 

(vi) The above choices of p and the scaling of ao is somewhat optimal 
in order to maintain the monotonicity. However, given G, they may not be 
optimal for the convergence of the scheme. In our numerical examples in 
Section 6 below, we may choose some different o"o and p. It is not clear to us 
how to choose p and ao so as to optimize the efficiency of the algorithm. 

Remark 3.7. This remark concerns the degeneracy of G. 

(i) Unlike [15], we do not require G 7 > \a^ and thus F can be degenerate. 
This is possible mainly because we use bounded trinomial tree instead of 
unbounded Brownian motion. 

(ii) Our condition D[G y ] > aid > is slightly weaker than the nonde- 
generacy requirement of [15]. In general, G 7 can be degenerate. 
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(iii) When a — 0, one can approximate the generator G by G e := G+ea^ : 
7 and numerically solve the corresponding solution u e . By the stability of 
viscosity solutions we see that u e converges to u locally uniformly. 

(iv) By imposing some additional assumptions on G z and by using the 
so called weak monotonicity, we may weaken Assumption 3.3 (ii) to A < Ag. 
We refer to [15] Assumption F (iii) for more details. 

(v) Motivated from pricing Asian options, in a recent work Tan [27] in- 
vestigated the numerical approximation for the following type of PDE with 
solution u(t, x, y)\ 

-d t u(t, x, y) - G(t, x, y, u, D x u, D 2 xx u) - H(t, x, y, u, D x u, D y u) = 0, 

where G is nondegenerate in D^ x u, but the PDE is always degenerate in 
D 2 yy u. 



3.3 Stability 

Given the monotonicity, one may prove stability following standard argu- 
ments. 

Lemma 3.8. Let Assumptions 2.1 and 3.3 hold. Then for any h e (0, ho], 
satisfies the stability in Theorem 2.2 (iii). 

Proof. First, it follows from Lemma 3.5 that satisfies the monotonicity. 
Denote C n := sup^d \u h (t n ,x)\ and d := sup (M)eMi+i)xKd \u h (t,x)\, i = 
n — 1, • • ■ ,0. Since g is bounded, we see that C n < C. We claim that 

Q < (1 + Ch)C i+l + Ch. (3.13) 

Then by the discrete Gronwall Inequality we see that 

max Ci < C(l + Ch) n [C n + nh\ < C. 

0<i<n— 1 

This proves the lemma. 

We now prove (3.13). Let t G [U,ti + i) and denote h :— ti — t < h. Similar 
to (3.11), one may easily get 



u h (t,x) 



E 



u h {t i+1 ,x + Vhao^Ii+i +hF(u,x, 0,0,0) (3.14) 



where, for some deterministic F y (ti), F z (ti), F 1 {ti) defined in an obvious way, 



H+l 



1 + h 



F y (U) + F g {U) ■ K^) + F 7 {ti) : K 2 (0 
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The monotonicity in Lemma 3.5 exactly means Jj + i > 0. Then, noting that 
F (ti, x, 0, 0,0) =G(ti,x, 0, 0, 0) is bounded, 



+ Ch < C i+1 E[I i+1 ] + Ch. 



\u h (t,x)\ < E \u h (t i+1 ,x + Vh,a ^)\I i+1 

By (3.6), we see that 

E u [I i+1 ] = l + hF y (t i ) <1 + Ch. 

Then 

\u h (t,x)\ < (l + Ch)C i+1 + Ch<(l + Ch)C i+1 + Ch. 
Since (t,x) is arbitrary, we obtain (3.13). □ 

3.4 Boundary condition 

It can be seen that the boundary condition holds if we can prove that 
\im (t/ y A) ^ iT:X:0) u h (t',x') = g(x) at T. 

Lemma 3.9. Let Assumptions 2.1 and 3.3 hold, then 

\u h (t k ,x) - g{x)\ < C(T-t fc )* 

Proof. Fix (tk,x). Let £ 3 , j = k + 1, • • • , n are independent (^-dimensional 
random variables such that each has independent components with distri- 
bution (3.2), and denote 



XI := x, XI := X%_ x + y/ha £ ?, J = k + ^ 



n. 



i=k+l 



Then it is clear that, denoting J 7 ™ := cr(£*, k + 1 < i < j), 
u^t^XU = VJu^,^)] +hF(t j _ u X?._ 1 ,E t ._ 1 [u h (t j ,X?.)] 

Similarly to the proof of Lemma 3.8, we have 

tifcfe-x.AgJ = E^^fe,^)/,] + / i F(t i _ 1 ,X£_ l , 0,0,0), 
where, by abusing the notation / slightly, 



1 + h 



F y {t^ x ) + F z {tj-x) ■ Me) + f^vx) : K 2 (e) 



> o, 
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and F y (tj-i), F z (tj-i), F y (tj-i) are defined in an obvious way. Denote 

J k := 1 and Jj := ir j=k+1 Ij, i = k + 1, • • • , n. 
Recalling Uh(t n ,x) = g(x), by induction we get 



n-l 



u h (t k , x) = u h (t k , XI) = E [g{Xl)J n + JjF(tj, X* , 0, 0, 0) 

j=k 

Since g is bounded and uniformly Lipschitz continuous, we may let g e be a 
standard smooth molifier of g such that 

||<7 e -<7||oo<Ce, II^IU < C and < Ce~\ (3.15) 

Then, noting that F(t, x, 0, 0, 0) is bounded, 

\u h (t k ,x) - g(x)\ 



n-l 



< 



\E[g £ (Xl)J n ] -g £ (x) 



CeE 



J 7) 



+ C7iE[^J i ] +Ce 



j=k 



n n 

= | E [9e(*Z)Ji ~ geiXZJJi-i] | + CeE[J n ] + ChE[j2 J 3 ] + ^ 



j=fc+l 



j=k 



i=fc+l 

n-l 



+CeE 



j=k 



Ce. 



Since E^_ 1 [/ 7 -] = 1 + hF y (tj-i), we have 

Eti-JA] -1 <Ch and E[JJ < (1 + C7i) l ~ fc < C. (3.16) 

Thus 

\u h (t k ,x) - g(x)\ (3.17) 

n 

< | £ E[[^(X«) - fc^J] Jj | + C(n - + Ce. 
i=fe+i 

Moreover, for some appropriate J-^-measurable X™, 
g e (X») - g e {XU = VhDg e {XU • a C + ^D 2 g e (X^ : (aoC^oCf ■ 
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Since is bounded, we have 



E u _, [D 2 g e {XD : (trofX^f/J < Ce" 1 !^ [/<] < C^ 1 



Then 



V, - < C/i[l + e" 1 ]. 

Plug this into (3.17) and recall (3.16), we have 

\u h (t k ,x) - g(x)\ < C(n- + e"" 1 ] +Ce. 



Note that (n — k)h = T — t k . Set e := y/T — t k , we obtain the result. □ 
3.5 Convergence Results 

First, combine Lemmas 3.2, 3.5, 3.8, and 3.9, it follows immediately from 
Theorem 2.2 that 

Theorem 3.10. Let Assumptions 2.1 and 3.3 hold. Then Uh converges to u 
locally uniformly as h — )■ 0. 

We next study the rate of Convergence. We first consider the case that 
u is smooth. Let C^ 4 '([0,T] x M. d ) denote the set of bounded functions u 
such that the following derivatives exist and are continuous and bounded: 

d t u, D x u, D 2 xx u, d 2 t u, d t D x u, d t D 2 xx u, D 3 xxx u, D A xxxx u. 

Theorem 3.11. Let Assumptions 2.1 and 3.3 hold and h G (0, ho). Assume 
further that u G C*j; 4 '([0,T] x M. d ), and F is locally uniformly Lipschitz con- 
tinuous in x, locally uniformly on (y,z,j). Then there exists a constant C, 
independent of h (or n), such that 

\uh{t,x) — u(t,x)\ < Ch for all (t,x). 

Proof. Again, since h < ho, it follows from Lemma 3.5 that satisfies 
the monotonicity. Denote C n := sup xgM d \uh(t n , x) — u(t n ,x)\ and C{ := 
su P(t,x)e[ti,t i+1 )xW \ u h(t,x) - u(t,x)\, i = n - 1, • • • ,0. We claim that 



d < (1 + Ch)C l+l + Ch 2 . 



(3.18) 
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Since C n = 0, then by the discrete Gronwall Inequality we see that 
max d < C(l + Ch) n \C n + nh 2 } < Ch. 

0<i<n-l 

This proves the theorem. 

We now prove (3.18). Similar to the proof of Lemma 3.8, we shall only 
estimate \uh(ti,x) — u(ti,x)\, and the estimate for the general \uh(t,x) — 
u(t,x)\ is similar. For this purpose, recall (3.4), (3.5) and define 

u h (ti,x) := [V°u(t i+1 ,-)](x) (3.19) 
+hF(u,x, [V°u(t i+1 , -)}(x), [V l u(t i+l , -)]{x), [V 2 u(t l+1 , -)](x)Y 



We note that the right side of above uses the true solution u, instead of Uh 
in (2.3). It is clear that 

\u h (U,x) - u(U,x)\ < \u h (ti,x) - u h (U,x)\ + \u h (ti,x) - u(U,x)\. (3.20) 

Compare (2.3) and (3.19), by the first equality of (3.12) we have 

u h (ti,x) - u h (ti,x) 
= E^Uh-ulit^x + VhaoOil + hlFy + Fz-K^O+F^ : K 2 (£)] 

Then it follows from similar arguments in the proof of Lemma 3.8 that 

\u h (U, x) - u h (U, x) | < (1 + Ch)C i+l . (3.21) 

Next, since u G C^([0,T] x ]R d ), applying Taylor expansion and by (3.3) we 
have 

[V°u{t i+1 , ■)] (x) = u{ti, x) + d t u(t h x)h + -D 2 u{U, x) : a 2 + 0(h 2 ); 

[V^iU+u •)] (x) = Du{ti,x) + 0{h); 
[V 2 u(t i+1 , •)] (x) = D 2 u(t h x) + 0(h); 



Then 

h 

u h (ti, x) - u(t u x) = d t u(ti,x)h + -D 2 u(ti,x) : al + 0(h 2 ) 



+hF^-, u + 0(h), Du + 0(h), D 2 u + O(h)^ (t u x) 
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Note that u satisfy the PDE (1.1) and recall (3.1), then 

iih(ti,x) — u(ti,x) =0(h 2 ) + hx 
F (-,u + 0(h), Du + 0(h), D 2 u + 0(h)) - F(-,u,Du,D 2 u) (t h x). 

Since u and its derivatives are bounded, and F is locally uniformly Lipschitz 
continuous in x, then we have 



u h (ti,x) - u(t u x) 



< Ch 2 . 



Plug this and (3.21) into (3.20), we obtain 

sup \u h (U,x) - u(U,x)\ < (1 + Ch)C i+ i + Ch 2 . 

x 

Similarly we may estimate sup^. \uh(t,x) — u(t,x)\ for t G (ti,t i+ i), and thus 
prove (3.18). □ 

We finally study the case when u is only a viscosity solution. Given the 
monotonicity, our arguments are almost identical to those of [15] Theorem 
3.10, which in turn relies on the works Krylov [18] and Barles and Jakobsen 
[2]. We thus present only the result and omit the proof. 

The result relies on the following additional assumption. 

Assumption 3.12. (i) The generator G is of the Hamilton- J ocobi-Bellman 
type: 

G(t, x, y, z, 7 ) = inf { \a a (a a ) T (t, x) : 7 + b a (t, x)y + c a (t, x) ■ z + f a (t, x) } , 

where the junctions a a , b a , c a and f a are uniformly bounded, and uniformly 
Lipschitz continuous in x and uniformly Holder-^ continuous in t, uniformly 
in a. 

(ii) For any 5 > 0, there exists a finite set {o;j} i= j such that for any 
a e A: 

inf (\a a - a a % + \b a - b a % + \c a - c a % + \f a - f ai U < 5. 

l<i<M 6 

We then have the following result analogous to [15] Theorem 3.10. 
Theorem 3.13. Let Assumptions 2.1 and 3.3 hold and h e (0, ho). 

(i) under Assumption 3.12 (i), we have u — Uh< Ch 1 ^ , 
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(ii) under the full Assumption 3.12, we have — C/i 1 / 10 < u — u h < Ch 1 ^ . 

Remark 3.14. The arguments in [15] relies heavily on the viscosity proper- 
ties of the PDE. Very recently Tan [26] provides a purely probabilistic argu- 
ments for HJB equations. His argument works for non-Markovian setting as 
well and thus provides a discretization for second order BSDEs. Moreover, 
under his conditions he shows that \v,h — u\ < Ch», which improves the left 
side rate in Theorem 3.13 (ii). 

4 Quasilinear PDE and Coupled FBSDEs 

In this section we focus on following G which is quasilinear in 7: 
G = \[< 7<7T \{ t ^ x ^y) '■ 7 + b(t,x,y,a(t,x,y)z) ■ z + f(t,x,y,a(t,x,y)z).(4 : .l) 

Here / is scalar, b is IR d - valued, and a is ~R dxm - valued for some m. In this 
case the PDE (1.1) is closely related to the following coupled FBSDE: 

X t = x+ [ b{s,X s ,Y s ,Z s )ds+ [ a(s,X s ,Y s )dW s ; 

Jo t Jo ,t (4.2) 



Y t = g(X T )+ [ f(s,X s ,Y s ,Z s )ds- f Z 
Jt Jt 



dW K 



Here W is a m-dimensional Brownian motion, and the solution triplet (X, Y, Z) 
takes values in M d , R, and M m , respectively. Due to the four step scheme of 
Ma, Protter, and Yong [19], when the PDE (1.1) has the classical solution, 
the following nonlinear Feynman-Kac formula holds: 

Y t = u(t,X t ), Z t = a(t,X u u(t,X t ))Du(t,X t ). (4.3) 

The feasible numerical method for high dimensional FBSDEs has been 
a challenging problem in the literature. There are very few papers on the 
subject, most of which are not feasible in high dimensional cases. To our 
best knowledge, the only work which reported a high dimensional numerical 
example is Bender and Zhang [7]. 

Our scheme works for quasilinear PDE as well, especially when aa T is 
diagonal (or close to diagonal in the sense that 9 is small), and thus is ap- 
propriate for numerically solving FBSDE (4.2). We remark that the do we 
will choose is different from o~(t,x,y) in (4.1), and the F defined by (3.1) is 
different from /. We shall present a 12-dimensional example, see Example 
6.4 below. 
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One technical point is that the G in (4.1) is not Lipschitz continuous in 
y, mainly due to the term |[cx<7 T ](t, x, y) : 7. This can be overcome when the 
PDE has a classical solution u. In general viscosity solution case, one may 
construct some approximating PDE which has classical solution and then use 
the stability of viscosity solutions, in the spirit of Remark 3.7 (hi). 

Theorem 4.1. Let G take the form (4.1). Assume 

(i) o~,b,f,g are bounded, continuous in all variables, and uniformly Lip- 
schitz continuous in (x,y,z); 

(ii) Assumption 3.3 holds; 

(111) The PDE (1.1) has a classical solution u e Cf ] ([0,T] x R d ) (as 
introduced in Theorem 3.11). 
Then \uh — u\ < Ch when h is small enough. 

Proof. We follow the proof of Theorem 3.11. Define Cj, i — 0, • • • , n and Uh 
as in Theorem 3.11 and again it suffices to prove (3.18). 
We first estimate \iih(ti,x) — Uh{U,x)\. Denote 

<f h (x) := u h (t i+ x,x), ip(x) := u(t i+1 ,x), ip:=(p h -ip. 

Then 

u h (U,x) - u h (ti,x) 
= V°iP(x) + hF(t h x, V^hix)^ 1 ^)^ 2 ^] 

-hF(t h x, V°(p(x),V 1 cp(x),V 2 ip(x / 



V°ip(x) - ^a 2 : V 2 ip(x) + hG(t h x,V ip h (x),V 1 ip h (x),V 2 i Ph (x] 



-hG[ti,x, V ^^),!) 1 ^^),!) 2 ^^ 
V°ip(x) + h \c l i)(x) + C 2 i)(x) + C 3 tp(x) 
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where 

£ 2 *fj{x) 
C 3 ip(x) 



--a 2 : V 2 ^{x) + ^[aa T }(t h x,V°<f h (x)) : V 2 ^{x) 
+b(U, x, V°ip h (x), a(ti, x, T>°tp h (x))T> 1 tp h (x)) ■ V^x); 

b(ti, x, V°(p h (x), a(U, x, V Q ip h (x))V l <p h (x)) 
-b(ti, x, V°<p(x), <r(ti, x, V°(p(x))V l (p(x))] ■ X?V(^) 



+ 



f(U, x, V ip h (x), a(ti, x, V Lp h (x))V 1 Lp h (x)) 
-f(U, x, V°ip(x), a(U, x, V°<p(x))V l <p(x)) 



Let 77 denote a generic function with appropriate dimension which is uni- 
formly bounded and may vary from line to line. Since u is smooth with 
bounded derivatives, one may easily check that T>°ip, V 1 ip, V 2 ip are bounded. 
Then 



C l ijj{x) 

£ 2 ^{x) 
C 3 *p(x) 



-[<7a T }(u,x,v tp h (x)) ■ v 2 ^{x) - ^ 2 : v 2 4>(x) + t) ■ x>V0*0; 

r)V°ijj(x) 
r}V if)(x) 



+r] ■ 



a(t i ,x,V°(p h (x))V 1 ip h (x) - a(t u x,V°ip(x))V 1 (p(x) 



r]V ijj{x) + r] ■ a(U, x, V°ip h (x))V 1 ^(x) 



a(ti, x, V°ip h (x)) - a(U, x, V°(p(x)) V 1 ^) 



Then 



= f]V°ip(x) +r] ■V 1 ip(x). 

u h (U,x) - Uh(U,x) 

V°ip(x) + h \r]V ijj{x) + 7] ■ V 1 tjj(x) 



+ 



E 



2 L 
h, 



: V 2 ip(x) 



^(x + Vha O [l + h V + h V - + 2 - ^ : ^(0] 



Now following the same arguments as in Lemma 3.5 we see that, for h small 
enough, 

1 + hT] + hT] . Kl ($ + ^[ aa T _ ^ . K ^ > 
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Then it follows from the arguments in Theorem 3.11 that 



\u h (ti,x) - u h (U,x)\ < (1 + Ch)C i+1 . 
Similarly we may prove 

\u h (ti,x) - u(U,x)\ < Ch 2 . 
Thus we prove (3.18) and hence the theorem. □ 

5 Monte Carlo Simulations 

In this section we discuss how to implement the scheme. Fix n > 1 and 
x G M. d , our goal is to numerically compute Uh(t ,x). For this purpose, we 
need to apply the scheme (3.4)-(3.5) on trinomial trees. Let £ l , % — 1, • • • , n be 
independent and each takes values in M. d and has independent components 
with distribution (3.2). Define the trinomial tree: 

X^:^x + VhJ2^e, i = 0,-..,n. (5.1) 
We next define Y£ := g(X^J, and for i — n — 1, • • • ,0, 

+hF(u, XI E ti [Y- + J, E tj Kr +1 K 1 (f +1 )], E ti [Y t lK 2 (C +1 )}) , 
where the kernels K\ and are defined in (3.5). Then it is clear that 

YZ = u h (ti,XZ). (5-3) 

In particular, Uh(to,x) = Y£. So it suffices to compute Y£. 

When the dimension is low, say d < 3, we may generate the whole tri- 
nomial tree, which consists of Yl^o^ + n °des. Then we may compute 
the exact value of Y n at each node, where the conditional expectation in 
(5.2) is computed by the weighted average. This method is very efficient and 
the result is deterministic. It is in fact comparable to the standard finite 
difference method. 

However, in above method the number of nodes grows exponentially in 
dimension d. When d is high, this method fails to work. As standard in 
the literature of BSDE numerics, in high dimensional case we use Monte 
Carlo approach to approximate the value of Y7\ In particular, we shall use 
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the least square regression method as introduced in Gobet, Lemor and Waxin 
[17]. We remark that Monte Carlo method is less sensitive to dimensions. For 
example, it can be seen in next section that we can use Monte Carlo method 
to approximate a 12- dimensional PDE with 160 time steps and 13333333 
paths, while for finite difference method 13333333 paths is insufficient for 2 
time steps when d = 12. 

Roughly speaking, for each t i( fix appropriate basis functions (f>j(ti,-), 
j = 1, • • ■ , J. We approximate conditional expectations E ti [Y^ i£jfe(£* +1 )]j 
k — 0,1, 2, in (5.2) by their least square regression on the linear span of 

{^j(^ x t!)}i<j<J- that is, by Ei=i%<M*i>^t") where { a j}i<j<J are a ( x u) 
measurable random variables minimizing 



j 

E[\J2^MU,XZ)-Y t lK k (C +1 



We next simulate L-paths of X n and use them to approximate the coefficients 
ctj. We refer to [17] for the details as well as the error analysis. 

We note that there are three types of errors in this algorithm: the time 
discretization error, the least square regression error, and the Monte Carlo 
simulation error. The first error is already analyzed in Sections 3 and 4. The 
third error is small when L is large, due to the Central Limit Theorem. The 
second error relies heavily on our choices of the basis functions. While there 
are some studies on how to choose good basis functions, see e.g. Bender and 
Steiner [6], overall speaking this is still an open problem. Usually g and its 
derivatives (when exist) are good candidates of basis functions. 

In our numerical examples in next section, we shall focus our attention to 
the discretization error and the simulation error. That is, for the examples 
where we know the true solution, we will choose basis functions whose linear 
span contains the true solution, and thus the regression error vanishes. We 
shall leave the choice of basis functions for future study. 



6 Numerical Examples 

In this section we apply our scheme to various examples. 

6.1 Examples under Monotonicity Condition 

In this subsection we consider examples with diagonal G 1 and we shall always 
choose do diagonal, so = and thus Ag = oo in Assumption 3.3. Therefore, 
we have no constraint on A. 
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We start with a 3- dimensional example for which we can compute its 
values over the trinomial tree by using the weighted averages. 

Example 6.1. A 3 -dimensional Fully Nonlinear PDE 



-d t u - ~ sup [{a 2 I d ) : D 2 u] + f{t,x,u,Du) = 0, on [0,T) 
u(T, x) = sin(T + x\ + ... + xj), on M. d , 

where < a < a are both in M., and 



x 



(6.1) 



f(t, X,y,z) = ^l>2 Z i- 2 a<f<a ' 
i=l ~ ~ 

It can be checked that (6.1) has a classical solution: 

u(t, x) = sin(t + xi + ... + Xd), (6.3) 

with which we can verify the convergence of our numerical approximation. 

To test its convergence under different nonlinearities, we assume that 
d = 3, a = 1, a = \/2, or Supposing that T = 0.5 and xq = (5, 6, 7), 
we know the true solution is u(0, xo) = sin(5 + 6 + 7) ~ —0.750987. 

According to our scheme, when G 7 is diagonal, 8 = 0, which implies that 

we can choose the following parameters: A = p = min j (a-ix^-i) ' l} ' 

a p — 2, a = 2p (A-i)+a ' anc ^ a ° = ^ e remar k that A = 2, 4, 6 respec- 

tively, which violates the contraint (2.8) of [15]. Denote the number of time 
partitions by n. By applying the weighted average method we can obtain 
the results in Figure 1, where the cost in time increases from 0.1 second to 
800 seconds exponentially as n increases from 20 to 160 linearly. The table 
in Figure 1 contains the numerical solutions when a 2 = 2 exclusively, while 
the graph depicts the errors under three different choices of o 2 . 

As we can see from Figure 1, the rate of convergence is approximately 
C ■ h, whereas the C depends on the structure of G. Therefore, our scheme 
works generally for large A when G is diagonal or diagonally dominant with 
a small 9. 

In Figure 2 we compare the convergence of our scheme with that of finite 
difference method by fixing q_ = 1, a = y/2. It can be seen that our result 
converges slightly slower than, but is comparable to, the finite difference 
method in solving low dimensional problems. 

To see more of our scheme in extreme condition, we assume a = 0. Then 
we truncate G 7 from below with a positive definite matrix eld > 0. That is, 
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n 


Ours 


F.D. 


20 


-0.72984 


-0.76420 


40 


-0.74028 


-0.75785 


60 


-0.74382 


-0.75562 


80 


-0.74667 


-0.75447 


100 


-0.74560 


-0.75379 


120 


-0.74738 


-0.75332 


140 


-0.74790 


-0.75300 


160 


-0.74829 


-0.75274 


Ans 


-0.75099 


-0.75099 



0.02 



2 0.015 



0.01 




Figure 2: Comparison with finite difference method in Example 6.1. 
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n 


Approx. 




a 2 = 2 


20 


-0.76285 


40 


-0.75705 


60 


-0.75508 


80 


-0.75401 


100 


-0.75339 


120 


-0.75297 


140 


-0.75269 


160 


-0.75247 


Ans. 


-0.750987 




Figure 3: Convergence of a degenerate PDE truncated in Example 6.1. 
we approximate (6.1) by the following nondegenerate PDE: 



-d t u - \ sup £<CT<7f 



[a 2 I d ) : D 2 u 



+ f(t,x,u,Du) = 0, < t < T, 



u(T, x) = sin(T + x x + ••• + x d ), on M d , e = 0.01, 

where / is given by (6.2). Figure 3 shows the feasibility of truncation in 
dealing with a — 0. □ 

Our main motivation is to provide an efficient algorithm for high dimen- 
sional PDEs. At below we test our scheme on a twelve dimensional example, 
for which we shall use the the regression-based Monte Carlo method sug- 
gested in [17]. 

Example 6.2. A 12- dimensional example 

Consider the PDE (6.1) with d = 12, a= 1, a = V2, 



f{t,x,y,z) 



d 



cos 



(t + \^ Xi) — - inf _ ( cr 2 sin(t + \^ Xi) 



1=1 



(6.4) 



i=i 



The true solution is u(t,x) = sm(t + Yli=i x i) a g am - As explained in 
Section 5, in this paper we want to focus on the discretization error and 
simulation error, so we rule out the regression error and test our algorithm 
by using the following perfect set of basis functions: 



x. 



sin(T + xi + ... + x d ), cos(T + x\ + ... + x d ) 
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n 


L 


K 


Avg(Ans.) 


Var(Avg.) 


cost (in seconds) 


2 


2083 


160 


0.659639 


3.53 x 10" 


-6 


4.48 x 10~ 2 


5 


13021 


64 


0.562635 


1.99 x 10' 


-6 


1.46 x lO" 1 


10 


52083 


32 


0.546598 


8.41 x 10" 


-7 


1.17 x 10° 


20 


208333 


16 


0.530432 


8.04 x 10- 


-7 


1.08 x 10 1 


40 


833333 


8 


0.521343 


2.25 x 10- 


-7 


9.11 x 10 1 


80 


3333333 


4 


0.519701 


1.21 x 10- 


-7 


7.28 x 10 2 


160 


13333333 


2 


0.517363 


6.17 x 10- 


-8 


5.86 x 10 3 


r 


?rue Solution 


0.513978 




n 



Figure 4: Numerical results of a 12-dimensional example in Example 6.2. 

To test the result, we fix T = 0.2 and Xq = {1, 2, 12}, which implies that 
the true solution is sin(78) = 0.513978. Assuming that we repeat K identical 
and independent tests, and we sample L paths in each test. Denoting the 
average of the results in K tests by Avg(Ans.), and the variance of this 
average by Var(Avg.), we can obtain the results in the table and figure in 
Figure 4, where we conducte fewer tests for larger L, because the results are 
stable enough to draw our conclusion. 

It can be seen from Figure 4 that the error shrinks slightly slower than 
O(h), which is due to the simulation error. Hence we want to explore the 
influence of simulation error by using all the parameters as above but fixing 
n = 40, K = 2, d = 12, T = 0.2, n = 40, a = 1, a = y/2. we increase L, the 
number of paths sampled, to see how the error reduces in Figure 5. While 
the variance and error decrease with more paths sampled, the cost in time 
increases linearly with respect to L from 8 seconds to 1400 seconds in Figure 
5. 
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L 


Approx. 


83333 


0.543643 


166667 


0.526979 


416667 


0.523897 


833333 


0.524683 


1666667 


0.521531 


333333 


0.521017 


6666667 


0.520083 


13333333 


0.518607 


Ans.: 


0.513978 



10 5 10 6 10 7 

L 

Figure 5: Relation between size of sample and errors in Example 6.2. 



We have seen that our scheme converges to the true classical solution if 
it exists. Meanwhile, if the PDE only has a unique viscosity solution, our 
scheme can render a converging result as well. 

Let / be zero in the Equation (6.1), then this equation has some un- 
known viscosity solution. However, our numerical results in Figure 6 still 
demonstrate a converging sequence. This can be also be observed from the 
decreasing differences between the numerical results. We shall remark though 
in this case our choice of basis functions may not be good. Again, we leave 
the analysis of the basis functions to future study. □ 

It is well known that Isaacs equations have a unique viscosity solution 
under mild technical conditions. We next test our scheme on the following 
Isaacs equation to see its performance. 

Example 6.3. A 12- dimensional Isaacs equation with unknown 
viscosity solution 



-u t -G(D 2 u) = 0, on [0,T) x 
u(T, ■) = sin(T + xi + ... + x d ), 



on 



where 



GH) : = > sup inf 

^0<«<lO<^<l 

i=i — — 

d 

= > inf sup 

^0<»<1 0<U <1 



a 2 (u,v)^u + f(u,v) 
a 2 {u,v)^a + f{u,v) 
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NR 5 - NR 2 


0.0334337 


NR 10 - NR 5 


0.0077685 


NR 20 - NR 10 


0.0076637 


NR m - NR 20 


0.0034146 


NR 80 - NR 40 


0.0012785 


NR im - NRso 
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Figure 6: Numerical results for a PDE with unknown viscosity solution in 
Example 6.2. 



and 



a 2 (u,v) := (1 + u + v), f(u,v):=- 




One can easily simplify G(7) as: G(j) = X/iLiS^T**) wnere 

{7«-J> 7«e(l,+oo), 
f + 7«e[-l,i], 
7« + |) 7« e -!)■ 

Therefore G(7) is neither concave nor convex when 7 = 0. Setting d=12, 
we assign arbitrary initial value xq = {xq }f =1 to inspect the outcome. One 
example tested here is Xq = 2m — T ~° ,5?r . 

Though the viscosity solution is unknown, our scheme still renders a 
converging numerical result in Figure 7. □ 

We next test our scheme for a 12-dimensional coupled FBSDE. 

Example 6.4. A 12-dimensional coupled FBSDE 
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Figure 7: A 12- dimensional Isaacs equation with unknown viscosity solution 
in Example 6.3. 

Consider FBSDE (4-2) with m = d = 12, a is diagonal, and 

1 1 d 

bi(t,x,y,z) := cos(y + Zi), au(t,x,y) := 1 + - sin (j,y~]xj + 

i=i 

f(t, x, y, z) := - sin(t + ^ x { ) [l + ^ sin {k^2 x i + v) 



=i i=i 

1 V-d 



1 + l sin (iJ2 j =i^j + y) 

d 

g(x) := sin(T + ^ 



— d cos(t + Xi) cos (y + cos(t + . 



i 

i=i i=i 

d 



Xi 
i=l 



The associated PDE (4.1) looks quite complicated, however, the coeffi- 
cients are constructed in a way so that u(t, x) := sin(£ + J2i=i x i) is the 
classical solution. Consequently, denoting X t := \^2,j =l Xl 



Y t = sin(t + dX t ), Z\ = cos(t + dX t ) 1 + - sin (X t + sin(t + dX t )) 

o 

satisfy the FBSDE. 

For PDE (4.1), we see that G 7 = |cr 2 is diagonal and | < a„ < | for 
each i. We note that / is not bounded and not Lipschitz continuous in y, 
however, since Z is bounded, then f(t,x,y,Z t ) is bounded and Lipschitz 
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Figure 8: A 12-dimensional coupled FBSDE in Example 6.4. 

continuous in y, and thus actually we may still apply Theorem 4.1. Set 
d = 12, T = 0.2, X = (2, 3, 4, 13), and apply the parameters specified 
before for our scheme. An approximation of Yq is shown in Figure 8, where 
the true solution Y t = sin(t + Ylt=i %t) nas value 0.893997 at t=0. □ 

6.2 Examples violating Monotonicity Condition 

In this subsection we apply our scheme to some examples which do not 
satisfy our monotonicity Assumption 3.3. So theoretically we do not know if 
our scheme converges or not. However, our numerical results show that the 
approximation still converges to the true solution. It will be very interesting 
to understand the scheme under these situations, and we shall leave it for 
future research. 

Example 6.5. A 12-dimensional PDE with a = 

Consider the same setting as Example 6.2 except that a — 0. 

Instead of truncating G 1 as we did at the end of Example 6.1, we will 
pick parameters p and o"o as if g_ were some small positive number. Then 
Assumption 3.3 is violated and our scheme is in fact not monotone. Never- 
theless, our numerical results show that our approximations still converge to 
the true solution, as presented in Figure 9. 
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Figure 9: A 12- dimensional example without monotonicity in Example 6.5. 



□ 

We next apply our scheme to the following HJB equation which is asso- 
ciated with a Markovian second order BSDEs, introduced by [9] and [25]: 

du I 

~ A7 - o SU P \-° 2 : ^ ~ /(*> X ) = °> 011 t°' T ) X fa k\ 

Ot 2 a<a<a ' (6.5) 

u(T,x)—g(x), on R d , 

When / = 0, this PDE induces exactly the G-expectation introduced by 
Peng [24]. We emphasize that, unlike in previous examples, here a, a, a G §> d 
are matrices and < a < a < a. In particular, G 7 is not diagonal anymore. 
We remark that one has a representation for the solution of this PDE in 
terms of stochastic control: 



u(0,x) = supE g{X%) + / f(t,X?)dt 



X? :=x+ a s dW ; 



where W is a rf-dimensional Brownian motion, and the control a is an W w - 
progressively measurable S d - valued process such that a < a < a. Due to this 
connection, these kind of PDEs and the related G-expectation and second 
order BSDEs are important in applications with diffusion control and/or 
volatility uncertainty. 

Example 6.6. A 10- dimensional HJB equation 

Consider the PDE (6.5) with g(x) = sin(T + X\ + + ^4) and appropriate 
f(t, x) so that 

u{t, x) = sm{t + x\ + — + ... + — ) 
is the true solution to the PDE. We set d = 10. 
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To begin our test, we select randomly an initial point X , and two 10- 
dimensional positive definite matrices a 2 and a 2 . Then we obtain the numer- 
ical results displayed in Figure 10, which shows that our method works well 
for high dimensional HJB equation. The parameters used in this example 



arc: 



X = (2.99, 3.05, 1.54, 1.89, 2.52, 1.10, 3.21, 1.64, 1.02, 1.80), 
which gives a true solution 0.75805, 



a 



I 1.18,-0.35,-0.29, 0.23,-0.52, 0.09,-0.09, 0.21, 0.25,-0.03 
-0.35, 2.84, 0.42,-0.23, 0.00,-0.03, 0.21,-0.38,-0.25, 0.73 
-0.29, 0.42, 1.54,-1.17, 0.32,-0.16,-0.64,-0.63,-0.35,-0.12 

0.23,-0.23,-1.17, 2.54,-0.30, 0.07, 0.30, 0.97, 0.43, 0.22 
-0.52, 0.00, 0.32,-0.30, 1.77, 0.25, 0.09,-0.39, 0.19, 0.13 
0.09,-0.03,-0.16, 0.07, 0.25, 2.13, 0.23, 0.82, 0.65, 0.42 
-0.09, 0.21,-0.64, 0.30, 0.09, 0.23, 1.79, 0.31, 0.06,-0.30 
0.21,-0.38,-0.63, 0.97,-0.39, 0.82, 0.31, 1.93, 0.14, 0.88 
0.25,-0.25,-0.35, 0.43, 0.19, 0.65, 0.06, 0.14, 1.39,-0.05 
-0.03, 0.73,-0.12, 0.22, 0.13, 0.42,-0.30, 0.88,-0.05, 1.76 



and 



a 



\ 



0.73, 
-0.21, 
-0.09, 

0.28, 
-0.18, 
-0.07, 
-0.07, 
-0.16, 

0.20, 
-0.22, 



-0.21,- 
1.63, 
0.15, 
-0.04,- 
-0.07, 
-0.30, 
-0.04, 
-0.40,- 
-0.19,- 
0.08,- 



-0.09, 
0.15,- 
1.06,- 
-0.80, 
0.18,- 
-0.28, 
-0.64, 
-0.66, 
-0.35, 
-0.06, 



0.28,- 
0.04,- 
0.80, 
1.31, - 
-0.07, 
0.57, 
0.19,- 
0.60,- 
0.69,- 
0.15, 



-0.18,- 
-0.07, 
0.18,- 
-0.07, 
0.38, 
0.10, 
-0.23, 
-0.04, 
-0.12, 
0.01, 



-0.07,- 
-0.30, 
-0.28,- 
0.57, 
0.10,- 
0.54,- 
-0.16, 
0.53, 
0.18, 
0.29,- 



-0.07, 
-0.04, 
-0.64, 
0.19, 
-0.23, 
-0.16, 
1.32, 
0.04, 
0.06, 
0.53, 



-0.16, 
-0.40, 
-0.66, 
0.60, 
-0.04, 
0.53, 
0.04, 
0.98, 
0.17, 
0.61,- 



0.20,-0.22 
-0.19, 0.08 
-0.35,-0.06 
0.69, 0.15 
-0.12, 0.01 
0.18, 0.29 
0.06,-0.53 
0.17, 0.61 
0.81,-0.11 
-0.11, 1.09 



One can check that a 2 > <r 2 because the smallest eigenvalue of a 2 — a 2 is 
0.0026, which is positive. □ 
Note that the PDE (6.5) involves the computation of svLp a<a<w [a 2 : 7]. 
We provide some discussion below. 

Remark 6.7. Let a 2 — cr 2 = LL T be the Cholesky Decomposition, namely L 
is a d x d lower triangular matrix. Then for any 7 e S d , we have 



sup [a 2 : 7] = a 2 : 7 + V 7^ 



a 2 <a 2 <W 2 



i=l 



where ^i, i = 1, ■ ■ ■ , d, are the eigenvalues of L Tr jL. 
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Figure 10: A 10- dimensional HJB equation in Example 
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Proof. Obviously, any a 2 G S d between ex 2 and a 2 can be expressed as a 2 = 
a 2 + A, where < A < LL T . Then < L~ X AL~ T < I d . We make the 
following eigenvalue decompositions: 

L~ x ALT T = UAU T , L T ^L = P^fP T , 

where UU T = PP T = 1^, and A and g diagonal matrices. It is clear that the 
diagonal terms of A are fij G [0, 1] and the diagonal terms of 7 are 7$. Denote 
Q := U T P. Then 

a 2 ; 7 -a 2 : 7 = A: 1= [ZT 1 AL~ T ] : [L T 7 L] = [£/if/ T ] : [L T 7 L] 
= i : [U T L T 1 LU] = A : [Q 7 Q T ] 

d d d d 

1=1 j=i 1=1 j=i 

Note that Xw=i ?ij = 1- Then by Jensen's inequality, 

„• : 7 - 2 * : 7 < E (E4V) + £ EE4¥ = Etf E«S = 

i=l j'=l i=l j=l j'=l i=l j'=l 

This proves the remark. 

Moreover, from the proof we see that the equality holds when 

hi = 1 {J2U^>°} and Q = Id - 

That is, U — P and thus a 2 = a 2 + LPAP T L T , where A is the diagonal 
matrix whose diagonal terms are <Zj = 1/A i>0 i. □ 

We remark that the above computation is in fact quite time consuming. 
At below we provide another example where G 1 is tridiagonal and the scheme 
is much more efficient. 

Example 6.8. A 10- dimensional example with tridiagonal struc- 
ture 

Consider the PDE (1.1) with 

G(t,x,y,z,i) := (sYLilH + S[i-j|=i v / 1+ ( 7ij)2 ) + f& x ) (6.6) 
g(x) := sin(T + xi + x 2 + ...x d ), 

and /(t, x) is a chosen function such that u(t, x) : = sin(£ + x% + ... + Xd) is 
the true solution of the PDE. 



33 



11 


L 


K 


Avg(Ans.) 


Var(Avg 


•) 


cost (in seconds) 


2 


2500 


160 


-1.47362 


1.0 x 10 


-5 


1.2 x W-' 2 


5 


15625 


64 


-1.15004 


1.7 x 10 


-6 


1.4 x lO" 1 


10 


62500 


32 


-1.06194 


9.1 x 10 


-6 


1.0 x 10° 


20 


250000 


16 


-1.04519 


2.1 x 10 


-6 


8.9 x 10° 


40 


1000000 


8 


-1.03326 


6.2 x 10 


-7 


7.1 x 10 1 


80 


4000000 


4 


-1.03092 


5.8 x 10 


-8 


5.9 x 10 2 


160 


16000000 


2 


-1.01910 


3.0 x 10 


-9 


1.4 x 10 4 




Figure 11: A 10-dimensional example with tridiagonal generator in Example 
6.8. 

In this case one may check straightforwardly that 

[G 7 ]a = 3 and [G 7 (t, x, y, z, 7)]^ = %1 — T , \i - j\ = I. 

When d — 10, this example is out of the scope of our monotonicity As- 
sumption 3.3. However, if we test it using T = 0.2, x = (1, 2, 10), the 
numerical results show that our scheme still converges to the true solution, 
sin(55) = —0.999755, as presented in Figure 11. 

We shall remark though that this example is computationally more ex- 
pensive than Example 6.2 because here we need to approximate 3d — 2 second 
derivatives. □ 
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